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Phonon-phonon interaction is systematically studied by nonequilibrium Green's function (NEGF) 
formulism in momentum space at finite temperatures. Within the quasi-particle approximation, 
phonon frequency shift and lifetime are obtained from the retarded self-energy. The lowest order 
NEGF provides the same phonon lifetime as Fermi's golden rule. Thermal conductance is predicted 
by the Landauer formula with a phenomenological transmission function. The main advantage of 
our method is that it covers both ballistic and diffusive limits and thermal conductance of different 
system sizes can be easily obtained once the mode-dependent phonon mean free path is calculated 
by NEGF. As an illustration, the method is applied to two one-dimensional atom chain models 
(the FPU-/3 model and the (^* model) with an additional harmonic on-site potential. The obtained 
thermal conductance is compared with that from a quasi-classical molecular dynamics method. The 
harmonic on-site potential is shown to remove the divergence of thermal conductivity in the FPU-/? 
model. 

PACS numbers: 05.60.Gg, 44.10.+i, 63.22.-m 
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I. INTRODUCTION 

More understanding of nanoscale thermal transport is 
required to solve heat dissipation problem, which be- 
comes important as the size of electronic device de- 
creasesi. In recent years, many experiments have 
been done to measure thermal conduction of nanostruc- 
tures^i^i^i^. At the nanoscale, the ballistic approxima- 
tion is usually a good starting point for thermal trans- 
port and a lot of theoretical research on ballistic ther- 
mal transport have been reported^iii^i^. However, the 
ballistic approximation would lead to unphysical results, 
such as infinite phonon mean free path, divergent ther- 
mal conductivity and zero temperature gradient. In most 
cases we need to go beyond ballistic limit to include 
effects of scattering for a more realistic consideration. 
In fact, phonon-phonon interaction is one of the sig- 
nificant factors for understanding and improving ther- 
mal transport properties. Since the sizes of nanostruc- 
tures are comparable to phonon mean free path, ther- 
mal transport in these systems is in the intermediate re- 
gion between ballistic and diffusive ranges. Nanoscale 
thermal transport has been studied through many ap- 
proachesiSiiiii^iiliiiii^ii&iilii^iaiSO. Nevertheless, no sat- 
isfactory method has been available to deal with thermal 
transport in the intermediate region. 

Our work aims to provide an efficient way to study 
ballistic-diffusive thermal transport including the effects 
of phonon-phonon interaction. Here we propose a new 
formalism of nonequilibrium Green's function (NEGF) 
to treat phonon-phonon interaction in momentum space 
at finite temperatures. It is known that NEGF is usu- 



ally used to study nonequilibrium problems, but its for- 
mulism can be applied to equilibrium systems as well. 
While Matsubara formalism^i is conventionally used for 
equilibrium Green's function calculation at finite temper- 
atures, NEGF is an alternative that also works for such 
situations. In our framework of NEGF, the analytical 
form of phonon self energy to any order can be easily 
derived and phonon frequency shift and lifetime are di- 
rectly related to the retarded self energy. The phonon 
frequency shift can also be accurately predicted by an 
effective phonon theory based on the ergodic hypoth- 
esis (equipartition theorem) ^^'^'^:^'* . For weak phonon- 
phonon interaction, however, frequency shift is usually 
very small and the major influence on thermal transport 
is from the finite phonon lifetime. It will be shown that 
NEGF of lowest order is equivalent to Fermi's golden rule 
when considering phonon lifetime. 

We also provide a new approach to study thermal 
transport using NEGF. Our NEGF formalism gives 
phonon lifetime and mean free path. Using a phenomeno- 
logical transmission function determined from mean free 
path in the Landauer formula, ballistic-diffusive thermal 
transport can be studie d^"'^^:^^ . Many works have been 
done on thermal transport by NEGF— liiiiSiii^i^S. In pre- 
vious treatments, the system is situated at the center as 
a junction with semi-infinite leads on the two sides. The 
central part can not be too large due to the constraint 
of computational capability. Here we compute phonon 
lifetime in a periodic system at equilibrium, and feed the 
lifetime information to phonon transmission function in 
the Landauer formula. Our approach is computation- 
ally more efficient, but with a less rigorous treatment of 



transmission function. 

In Sec. II, the general theory of NEGF is developed for 
phonon-phonon interaction and ballistic-diffusive ther- 
mal transport. As an application we employ the method 
to study two explicit models in Sec. III. A summary is 
made in Sec. IV. 



II. GENERAL THEORY 



two schemes are provided to solve Green's function. One 
is to use the equation of motion and recursive expan- 
sion rules. The other is to apply Feynman rules for self 
energy and solve the Dyson equation. Within the quasi- 
particle approximation, the retarded self energy directly 
corresponds to phonon lifetime. Next we will compare 
the NEGF method with Fermi's golden rule. Finally, we 
discuss how to study ballistic-diffusive thermal transport 
with the information of phonon mean free path. 



In the harmonic approximation, a crystal can be de- 
scribed in terms of non-interacting phonons. The concept 
of "phonon" remains valid when the anharmonic contri- 
bution is small compared with the harmonic. In this case, 
the quasi-particlc approximation can be made, and then 
the anharmonic effects give a complex shift to phonon 
frequency^^: the real part shifts the value of frequency 
while the imaginary part corresponds to phonon lifetime. 

One of the earliest method o^^'^^i'^" of calculating 
phonon lifetime is based on Fermi's golden rule and the 
Boltzmann-Peierls equation. However, those approaches 
can not be systematically improved. In contrast, the 
Green's function method provides a systematic way to 
consider phonon-phonon interaction at finite tempera- 
tures. Phonon-phonon interaction has been studied us- 
ing equilibrium Green's function metho d^'^^'^-'^i'^^'^'^i'^^'^^ . 
These works can be roughly divided into two groups ac- 
cording to the decoupling schemes used. For the first 
group, the temperature dependent Green's function with 
Matsubara representation is employed based on a gen- 
eral Wick's theorem which holds when the system is large 
enougb22i^i22i^. A fictitious imaginary time is used and 
an analytical continuation is needed to get the retarded 
Green's function^. For the other group, a decoupling 
scheme proposed by Valle and Procacci is used^i^^. In 
detail, in the equation of motion for the real time Green's 
function, thermal averages are replaced with those appro- 
priate for the harmonic Hamiltonian. The analytical ex- 
pression of phonon self energy including high order terms 
can be easily evaluated by a computer-aided technique, 
despite of the lack of a clear theoretical foundation. 

We will employ the NEGF method as described in 
Ref. [33. Its formalism has been generalized to treat arbi- 
trary nonequilibrium systems with arbitrary initial den- 
sity matrices by Wagner—. If interactions are adiabati- 
cally switched on and an initial non-interacting Hamilto- 
nian is used, Wick's theorem is still valid for the contour 
ordered Green's function^^. The adiabatic switch-on as- 
sumption leads to a loss of information related to initial 
correlation, but is still reasonable in most physical situa- 
tions. We will use this assumption and take the contour 
ordered Green's function to study phonon-phonon inter- 
action. 

In this section, we develop the general theory of NEGF 
to deal with phonon-phonon interaction and ballistic- 
diffusive thermal transport. At first, a general Hamil- 
tonian of anharmonic system is given, and the contour 
ordered Green's function of phonon is defined. Then, 



A. The Hamiltonian 

In a Taylor expansion of the potential energy $ at equi- 
librium configuration, terms higher than the second order 
constitute anharmonic Hamiltonian: 



Ha 



^ — Y^ $»l»2-i„"ilW<2 ■■•"»„■ (1) 
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n— 3,4,--- Ji ,J2,"" ;^Ti 



Herein, the system has N unit cells and each cell con- 
tains r atoms. We use i = {rii, i) (ui — 1,2,- ■■ , N) and 
i = {Ki,ai) {Ki = 1,2,--- ,r; ai ~ x,y,z) for conve- 
nience. Ui — -y/MsiXi, Mk- is the atom mass and Xi is 
the displacement. We expand Ui as 

"' = E^-7^4'(q)^«' 9 = (q'J9)' (2) 

where j^ = 1, 2, • • • , 3r denotes the phonon branch, and 
e^«(q) is the normal mode eigenvector. The phonon op- 
erator Aq can be expressed by the usual phonon annihi- 
lation operator a^ and creation operator aj. Introducing 
9= (-Qiig), we have 



"J ^ \\ 9 — ^^'^ ^ "«''' 1 ~ ■" 



2uj, 



(3) 



where tOq is the frequency of the normal mode q. Then 
the harmonic Hamiltonian has the form 



^0 = 5Z 9^9^9 + Y. \^W^-^^ 



(4) 



Aq = -%\ 



iy^iaq-al). (5) 

The anharmonic Hamiltonian becomes 

^^= J2 - J2 Pqiq2-qr.AgiA.---Ar., (6) 

n=3,4,--- qi,q2,--- i9re 

where the n-leg vertex is 

A(qi + q2 -I h q„) 

{n 
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x6-(qi)C(q2)---ff"(q«)- 



(7) 



A(q) = 1 if q is zero or a reciprocal lattice vector, oth- 
erwise A(q) = 0. 



B. Nonequilibrium Green's function 

We define the contour ordered Green's function of 
phonon as 



G. 



,Ar,T)^--{TrA,{r)A,,ir')), 



(8) 



where t and r' are defined in the complex plane, (• ■ ■ ) 
denotes the ordinary Gibbs statistical average, and TJ- is 
the contour-ordering operator. The contour runs slightly 
above the real axis from — oo to +00, and back to —00 
slightly below the real axis. From our definition, we have 



Gqq' (r, r') = 6g^q'Ggg{T, r'). 



(9) 



Note that only those Green's functions of the form 
Gqq{T,T') are nonzero because of the orthogonahty be- 
tween the normal modes. 

The mapping of the contour ordered Green's function 
onto four different normal Green's functions has been dis- 
cussedi^i^i22. A label u = ±1 is needed to distinguish 
whether r is on the upper branch or on the lower branch, 

G^g^,' {t, t') = lim^ Gqq, [t + ^ea, t' + tea'). (10) 

G++ ^ G\ G— = G\ G+- = G< and G-+ = G> . 
There are another two types of Green's functions: G"' 
and G" . The relations among these six Green's functions 
in time and frequency domains have been presented in 
Ref. M. 



C. The equation of motion 

Simple recursive expansion rules for the contour or- 
dered Green's function of phonon in coordinate space is 
originally proposed in Ref. |20. Differently, the phonon 
operators Aq used here are not Hermitian. However, we 
can show that the same recursive expansion rules can be 
used by our phonon Green's function defined in momen- 
tum space. 

We define a general n-point Green's function 

Gqiq2---q,^{Tl,T2, ' ' ' , T„) 
= - '-{%Aq,(T^)Ag,(T2) ■ • • A„(t„)). (11) 

The phonon operators satisfy 

Aq = -]:[Aq,H] =~wlAq 



/ J / , ^qq2q3---qn^q2^q3 ' ' ' ^Qn- i-'-^J 

n=3,4,--- q2,q3,--- ,qn 

From the definition, 

Fqlq2■■■q,^{'^l^T'2:■ ' ' , ^n) 
= Fqiq2---qr.^Tx,a25{tl - ^2) ' ' ' 5a^^a„5{ti - i„)CT""^ (13) 



The equation for the contour ordered Green's function is 

Q:^Gqq'{T,T') = Sq.q, S{t - t') - WgGq^q, {t, t') 



dr-)--- dTn 



- E E 

n=3,4,--- 132, ■■■ ^q 
X ■P'«92---9„('^:^2,--- ,Tn)Gq^...q^q,{T2,--- ,Tn,T'). (14) 

The equation for the unperturbed Green's function is 
92 



9r'2 



G'qq,{T,T') = ~Sq^q, Sir - t') - W^, G^, (t, t') . (15) 



Combining Eqs. P^ and fTS]) . the two-point Green's 
function can be expressed in terms of the free and higher 
order ones as 

Gqq>iT,T')^G"gg,ir,r')+ J2 E 

n=3,4,--- qi;q2;--- ,<2n 
dTldT2 ■■ -dTriG^ (t, Tl) 



Fa- 



q„{Tl,T2, 



' '''■n)Gq2---q„q'\T2, 



,r') . (16) 



Repeating the above procedures, we can get the equation 
for higher order Green's functions. 



G. 



91 92 ■■■9. 



{tiT2 



= iflGg^g^{Tl,T2)Gq^q^...q^{T3,T4, ■ ■ ■ ,r„) 

+ ihGg^g^{Ti,T3)Gq2q^...q^{T2,T4 ■ ■ ■ ,r„) -) 

+ iflGg^g^{Tl,Tn)Gq2q3...q„_l{T2,T3 ■ ■ ■ ,T„_l) 

+ E E [J /■■■Jdr[dr^---drL 

m=3,4,--- q[,q'2,--- ,9j„ 
X ^qiq[i'''ii'''l)Fq[q'^---q'^^{Tl,T2, ■ ■ ■ , r„) 



xG, 



92---9m92 



■9.(^2 



27 ' ' ' ; 'm 



m ; '''2 7 ' ' ' i^nj 



(17) 



Eqs. pop and p7|) arc the equations of motion for the 
contour ordered Green's function, which have the same 
form as those in Ref. l20|. So the same recursive expan- 
sion rules can be used here. One may conveniently im- 
plement these rules in a computer program and expands 
the Green's function to any order as one wish. 



D. Feynman rules for self energy 

The Dyson equation for our Green's function is 
Gqq{T,r') = G''ggir,r') 

+ I I dT4dT2G%{T,Tl)Y.qq{T4,T2)Gqq(T2,T'). (18) 

For a steady or equilibrium system, it is more convenient 
to treat problems in the frequency domain, 

+ ^ aic72G°--n^)S^r^(^)G^I"'(^)- (19) 
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FIG. 1: Feynman diagrams for phonon self energy. 



Feynman diagrams for phonon self energy are the same 
as those in Ref. 20;. Fig. [T] shows some lowest order Feyn- 
man diagrams. In the following we summarize the Feyn- 
man rules which are used to write the proper self energy 

(1) Draw all topologically distinct Feynman diagrams 
for the proper self energy with two terminals on the left 
and right separately. These diagrams should be con- 
nected and can not be separated into two parts by cutting 
a single line. 

(2) Draw an arrow on each line from left to right and 
label it with new variable (qi,uji). The line represents 

an unperturbed Green's function Gg.q. ' (wi). cFj,aji — 
±1 are the variables of the vertices the line connects, 
which will be discussed below. The variables of the two 
terminals are both {q,Lo). Using the relations between 

the six Green's functions, Gq.q\ ^' [ui) can be obtained 
from 



o — *'U"'" 



■i5f 



w„ 



(20) 



(3) Label the vertices: each internal vertex with new 
variable ai {cii = ±1), the two terminal vertices with 
a on the left and a' on the right. If only one terminal 
vertex exists, the self energy is zero when a ^ a' . Each 



n-leg vertex contributes a factor F^ 



91/91. 92/92, ■■■ ,9n/9i 



. Use 



qi in F if the line (qi,LUi) goes into the vertex and use qi 
in F if the line goes out of the vertex. The momentum 
conservation is automatically included in the factor. 

(4) At each vertex a factor of '2Tr6{J2i ^Un ~ 
Si ^iout ) is associated with energy conservation. The 
line {qii^ , Wi .„ ) enters the vertex and the line (gi„„t , i^i^^t ) 
leaves the vertex. The two terminal vertices only con- 
tribute one energy conservation factor. 

(5) Multiply all the internal a variables and a coeffi- 



cient 



1+ Y. N^ 2+ E 



c=(-l) ->-- (z) 

xn[("-i)r"i 



-N„ 



{hY>-^ 



-Af„ 



n>3 



S' 



(21) 



where Nn is the number of n-leg vertices and S is the 
symmetry factor of the Feynman diagram. 

(6) The sums or integrations are performed over all the 
internal variables. 

For example, the self energy of Fig. [TJc) is 

91---95 <T1 '' 

X 2t:5{uj — uJi — UJ2 — uj3)2ttS{uji + 102 — 1^4 — ^5) 



X -f'ggi 5293^91 929495^9394959" 

X 0°;-^-^ (..3X'^-r'(^4)G°;^-r'(^5). 



qiqi l^lJ<^g292 \^'^) 



(22) 



In equilibrium system, only one of the six kinds of self 
energies is independent. They have the same relations as 
the Green's functions. For example, if T,'^^{uj) — S| (w) 
is calculated, the retarded self energy can be expressed 
as 



^Ju;) = Re[^Ju;)]+^ 



Imp* (0 



where /(w) is the Bose-Einstein distribution. 



(23) 



E. The retarded self energy and phonon lifetime 

In the quasi-particle approximation, the phonon fre- 
quency tOq suffers a complex shift Aq — iTq due to phonon- 
phonon interaction, where F^ is the reciprocal phonon 
lifetime and 2Tq is the full width at half-maximum 
(FWHM) of the phonon peak in spectra. Then the fre- 
quency of the quasi-phonon is cuq = tOq + Aq. The re- 
tarded phonon Green's function can be written as 

1 



GUoj) 



(24) 



^2 - (^LUq + Aq - iVq)^ 

The Dyson equation for the retarded Green's function is 



GUto) 



1 



^^-U^^-^qico)- 

The quasi-particle approximation is valid if 

\Aq -iVql < UJq. 



(25) 



(26) 



With this condition, from Eqs. ([241) and ((25)) we get the 
relations 

Re[EJ-,(^,)] - 2ujqAq, (27) 

Im[l];-^(c.,)] - -2c.,r, = -^. (28) 



The imaginary part of the retarded self energy gives 
phonon hfetime. The condition of the quasi-particle 
approximation becomes |S^^(a;g)| ^ tj^ . Note that 
S''(-w) = [S''(cj)]* impUes Re[S''(-w)] = Re[I]''(w)] and 
Im[I]''(— w)] = — Im[I]'''(w)]. The retarded self energy of 
"w-independent" diagrams is real. To consider phonon 
lifetime, we only need to calculate "w-dependent" dia- 
grams. 



F. The NEGF method and Fermi's golden rule 

Fermi's golden rule is widely used to calculate phonon 
lifetime in previous studies. What is the relation and dif- 
ference between the NEGF method and Fermi's golden 
rule? This is the question to be answered in this subsec- 
tion. 

Let's consider three-phonon interaction at first. The 
lowest order self energy contributed by the three-phonon 
interaction is described by the Feymann diagram shown 
in Fig.[lja). Use the above Feymann rules , we can write 
down the corresponding "lesser" self energy as 



'Z192 •' •' 

The free "lesser" Green's function is 



Wl - LO2) 



(29) 



LUq 

(30) 



The imaginary part of retarded self energy can be ob- 
tained through the relation 



Imp< (^)] = 2f{u;)lu,[^Ju;)]. 



(31) 



Using Eq. ((28| . the reciprocal phonon lifetime can be 
expressed as 



r.-E 



Trft 



<7l92 



-iujqUJq^UJq^ 



|J^g.igJ'{[/KJ + /KJ + l] 



X [5{UJ - UJq^ -UJq^)~ 5{UJ + UJq^ + UJq^)] 

+ [/Kl)-/K2)]['5('^+^9l - ^q^) - ^{^ - '^qi +^92)]}- 

(32) 

On the other hand, we can use Fermi's golden rule to 
solve this problem and obtain the same result of recipro- 
cal phonon lifetime, which is presented in Ref.|40|. Fermi's 
golden rule is also equivalent to the lowest order NEGF 
for more than three-phonon interaction. This conclusion 
can be proved in a similar procedure as above and it will 
not be elaborated here. 

In summary, NEGF of lowest order is equivalent to 
Fermi's golden rule but the higher order terms of NEGF 
give additional corrections to it. The NEGF method pro- 
vides a more comprehensive and systematic way than 
Fermi's golden rule to treat phonon-phonon interaction. 



G. Thermal transport from phonon mean free path 

Thermal conductance of a system at temperature T 
can be given by the Landaucr formula as 



9(1) J >0) 



hiOqVq 



dficOq) ^ 



dT 



"91 



(33) 



where L is the length of system in the direction of heat 
flow {x direction here) , v^ is the phonon velocity in the x 
direction, fiujq) is the Bose-Einstein distribution, and S^ 
is the transmission function of phonon mode q = (q, J9). 
Eq. (I33p is valid from one-dimensional (ID) to three- 
dimensional (3D). Thermal conductivity has the form of 
K = (tL/S in 3D systems, where S is the cross section 
area. And for quasi-lD systems, we define k = uL. 

A phcnomenological formula is proposed to calcu- 
late phonon transmission function from the mean free 



s, = (i + L/g-i 



(34) 



with L = 



'q'q- 



When L <C /g, ^9 = 1, corresponding 



to the ballistic limit; when L ':^ Iq, ^q = Iq/L, thermal 
conductivity is 



E 

9(uJ>0) 



huj, 



1 5/K) ., 



'V dT 



q''q^ 



(35) 



where F = L in ID case and V — LS for 3D case. It re- 
produces the well-known Debye-Peierls formula for ther- 
mal transport. The formula covers the range from ballis- 
tic regime to diffusive regime^'^^. 



III. APPLICATIONS 

We apply the general theory developed to study ex- 
plicit models in this section. Firstly we introduce two ID 
atom chain models and investigate them by NEGF. Then 
a quasi-classical molecular dynamics (QMD) method is 
applied to the same models and thermal conductances 
from both methods are compared. Finally calculated re- 
sults and related discussion are provided. 



A. ID atom chain models 

The ID atom chain models have been intensively stud- 
ied. One of remarkable problems is in what condition the 
thermal transport of a ID system obeys Fourier law. The 
Fcrmi-Pasta-Ulam (FPU) /? model and the (p* model are 
two example models which are extensively used due to 
their simplicity. It is found that the thermal conduc- 
tivity in the FPU-beta model diverges with system size 
due to the momentum conservation, while it is conver- 
gent in the 0* model^'— — — . Previous study indicates 



that the external potential plays a determinant role for 
normal thermal conduction^i^. 

We will study the FPU-/3 model and the (p^ model 
with an additional harmonic on-site potential. In the 
following, these two models are denoted as model I and 
model II, respectively. The Hamiltonian of model I is 
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{Ui 



H+lj 






The Hamiltonian of model II has the form 
ill K , so Kq r. 



H 



E 



^ + y(wi-w*+i) +— "I 



M 4 
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They share the same harmonic Hamiltonian: 



Hn 



E 






"„,2 



(36) 

■ (37) 

(38) 



The dispersion relation of unperturbed phonon can be 
expressed as ojq ~ y^Ko + 2K{1 — cos qa). We set the 
lattice constant a — 1 A, K — 1.0 eV/(amu A^) and 
Kf) ~ O.li^. These parameters are chosen to be the same 
as those in Ref. [iJ: where model II is studied by the 
QMD method. The group velocity is defined as Vg = 
dojq/dq. Since a harmonic on-site potential is included 
in ifo, Vq=Q = and ujq ranges from 0.31 x 10^''s~^ to 
1.99 X IQi^s-i. 

For model I, we choose 13—1.0 eV/(amu^ A'') and the 
NEGF method provides converged result till T = 600 K. 
For model II, ^ = 1.0 eV/(amu^ A**) is also tried , but 
the NEGF method fails even for T = 100 K. To get a 
converged result, we select /z = 0.05 eV/(aniu^ A^) and 
the NEGF method would still work at the temperature 
higher than 1000 K. This is partially because a small an- 
harmonic parameter /i is chosen. The calculations indi- 
cate that the NEGF method which is based on perturba- 
tion expansion is only suitable for weak phonon-phonon 
interactions. 

The two models share the same set of Feynman dia- 
grams for self energy. Since their anharmonic parts are 
both quartic, only those vertices with four legs will ap- 
pear in the Feynman diagrams. The 4- leg vertex of model 
I has the form of 

2/3, 



93 +94)-;t[1 -cosgza (39) 



^gig2?394 = ^(91 + <?2 
— cos (73a — cos (74a -I- cos ((72 + 93)0 + cos ((72 + 94)0 
-I- COS ((73 -f qi)a - cos (92 + 93 + 94)0]. 



As qi (i = 1,2,3,4) approaches zero, Fq^q^q^q^ goes to 0. 
The 4-leg vertex of model II is 



9192 9394 



A(gi + 92 



93+94)^. 



(40) 



Only "w-dependent" diagrams of lowest orders are con- 
sidered in our calculation. These include one diagram of 
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FIG. 2: (Color online) Comparison of thermal conductance 
((t) from different methods: NEGF (red solid symbols) and 
QMD (black open symbols), (a) Thermal conductance of 
model I with the system size A'^ = 64, 256, 1024. Ballistic 
thermal conductance is also provided by NEGF. (b) Thermal 
conductance of model II with the system size A^ = 64, 1024. 



0(A'*) (Fig.^b)) and two diagrams of 0{\^) (Fig. [^c) 
and Fig. [Hd)). Using the Feynman rules developed 
above, the analytical form of self energy for these three 
diagrams can be derived, and then be used in numer- 
ical calculation. As described in the previous section, 
phonon lifetime and mean free path can be obtained from 
the imaginary part of the retarded self energy. Thermal 
transport properties can be given by the Landauer for- 
mula. 



B. Comparison of thermal conductance from 
NEGF and QMD 

Based on a generalized Langevin dynamics, the QMD 
method is developed to study quantum thermal transport 
in Ref. Jjj, where the considered system consists of a cen- 
tral junction part and two leads serving as heat baths. 
This method uses quantum heat baths derived from Bose- 
Einstein statistics and treats the central part classically. 
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FIG. 3: (Color online) Results of Model I. (a) Phonon lifetime Tq a,t T — 100 to 600 K. The nearly infinite lifetimes 
of very long wavelength phonons are not presented in the figure. (b) Thermal conductivity k with the system size 
N = W\ 10^ 10^ lO"*, 10^ 10^ lO'^, 10^° at r = 100 to 600 K. 
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FIG. 4: (Color onhne) Results of Model II. (a) Phonon lifetime r, at T = 100 to 600 K. (b) Thermal conductivity k with the 
system size iV = 10\ 10^ 10^ 10^, 10^ 10^° at T = 100 to 1000 K. 



In Ref. [ij, it is proved that the method can produce 
correct results both in quantum ballistic and classical 
diffusive limits. 

The QMD method is employed to give thermal con- 
ductances of the two models. 3 x 10* MD steps are used 
in the calculation with the time step of 10~^^s. The 
comparison of thermal conductance with that from the 
NEGF method is presented in Fig. [2 for models I and 
II. Both methods give essentially the same results at low 
temperatures. Deviations appear at high temperatures. 
The NEGF method gives larger thermal conductance for 
model I and smaller thermal conductance for model II. It 
is not unexpected that they give quantitatively different 
results. A quasi-classical approximation is made in the 
QMD method. The NEGF method makes perturbation 
expansion and neglects higher order terms which become 
important at high temperatures. Both methods provide 
only approximate results, and it is still unclear which 
method is superior. More work is needed to understand 
these differences. 



C. Results and discussion 

The results of model I and model II are pre- 
sented in Fig. [3] and Fig. 3] respectively. The mode- 
dependent phonon lifetime at different temperatures and 
the temperature-dependent thermal conductivity with 
different system sizes are shown . 

The two models with different types of anharmonic 
potentials give distinct properties of phonon lifetime. 
For model I, the translational invariant quartic poten- 
tial leads to zero retarded self energy and infinite phonon 
lifetime for infinite long wavelength mode, which is not 
shown in Fig. [3{a). Long wavelength phonon mode has 
long lifetime. In contrast, the long wavelength phonon of 
model II, which experiences large scattering due to the 
existence of quartic on-site potential, has short lifetime. 

Increasing temperature has two effects on thermal con- 
ductivity: shorter phonon mean free path and more 
phonon excitation. The first effect decreases and the 
second effect increases thermal conductivity. When the 
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FIG. 5: (Color online) Thermal conductivity k of different 
system sizes for model I and model II at T = 100, 300 K. 



system length is much smaller than the phonon mean 
free path, the first effect has minor influence and ther- 
mal conductivity would increase with increasing temper- 
ature. At high temperatures, for a system whose length 
is comparable to the phonon mean free path, the exci- 
tation of phonon modes becomes less important to ther- 
mal transport and then thermal conductivity would de- 
crease with increasing temperature due to the shortening 
of phonon mean free path. These are verified in Fig. [Hl^b) 
and Fig.UJb). For very large system at low temperatures, 
the thermal conductivity of model I decreases as the tem- 
perature increases, which is different from that of model 
II. This may be explained by the fact that the tempera- 
ture has larger infiuence on the phonon mean free path 
in model I than in model II at low temperatures. 

One of the advantages of our method is that once 
the mode-dependent phonon mean free path is obtained, 
thermal conductivity of different system sizes can be ob- 
tained easily. It would be interesting to show the size- 
dependent behavior of thermal conductivity, which is 
presented in Fig. [Sj When the system length is much 
smaller than the phonon mean free path (i.e., the bal- 
listic regime) thermal conductance does not change with 
system size, and thermal conductivity increases linearly 
with the system size. When the system length increases 
to much larger than the phonon mean free path (of or- 
der of 10*^ A in model II at 100 K) (i.e., the diffusive 
regime), thermal conductivity will cease to change with 
system size. A ballistic-diffusive transition with increas- 
ing system size is shown in Fig. [5] clearly. 

The FPU-/3 model, model I with Kq = 0, has been 
widely known to have divergent thermal conductivity. 
With an additional harmonic on-site potential, we get 
a finite bulk-limit thermal conductivity for model I as 
shown in Fig. [5] How does a nonzero Kq^ the harmonic 
on-site potential, induce different behavior of thermal 
conductivity? An analysis on long wavelength phonon 
modes will answer this question. As q approaches zero. 



ductivity contributed by long wavelength phonon modes 
would diverge with increasing system size. When Kq ^ 0, 
LOq -^ constant, Vq (x q, Tq (x q~^ and Iq oc q~^. The 
quadratic on-site potential breaks translational invari- 
ance and the phonon modes with very long wavelength 
have nearly zero group velocity. Though those phonon 
modes has nearly infinite mean free path, it can be shown 
analytically that the thermal conductivity contributed 
by them does not diverge with increasing system size for 
Kq =/= 0. Other phonon modes which have finite mean 
free paths also contribute a finite thermal conductivity. 
So finite thermal conductivity is obtained in bulk limit 
for model I. 

Comparing the results of the two models may give some 
information on how an anharmonic on-site potential in- 
fiuences thermal transport. A quartic on-site potential 
is included in model II. The anharmonic parameter fi of 
model II is 20 times as small as that of model I. How- 
ever, model II gives much shorter phonon lifetime for low 
frequency modes, which are very important for thermal 
transport especially at low temperatures. As shown in 
Fig. [51 the bulk-limit thermal conductivity of model II is 
much smaller than that of model I at 100 K. The compar- 
ison indicates that even a small quartic on-site potential 
can largely decrease the phonon lifetimes of long wave- 
length modes and significantly affect thermal transport 
especially at low temperatures. Our results confirm that 
the external potential plays a determinant role in thermal 
transport. 



IV. SUMMARY 

We have provided a new approach to study phonon- 
phonon interaction and ballistic-diffusive thermal trans- 
port by the NEGF method and the Landauer formula. 
A new formalism of NEGF has been developed to sys- 
tematically study phonon-phonon interaction in momen- 
tum space at finite temperatures in equilibrium. Using a 
phenomenological transmission function, which can be 
obtained from the mode-dependent phonon mean free 
path given by the NEGF formalism, the Landauer for- 
mula predicts thermal transport properties from ballistic 
region to diffusive region. Our approach is efficient for in- 
vestigating ballistic-diffusive thermal transport in weak 
interaction situations, where little additional computa- 
tional effort is needed when the system size changes. As 
an application, we have investigated two ID atom chain 
models. The results obtained are qualitatively agree with 
those by the QMD method. It is found that an additional 
harmonic on-site potential in the FPU-/3 model could re- 
move the divergence of thermal conductivity and a small 
quartic on-site potential can largely reduce the phonon 
lifetimes of long wavelength modes. The results confirm 
that the external potential plays an important role in 
thermal transpor t—'^ . 
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